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METHOD AND SYSTEM FOR ENHANCING SOLUTIONS TO A SYSTEM OF 

LINEAR EQUATIONS 

5 

This invention was made with U.S. Government support under NIH grant number ROl- 
CA66184. The U.S. Government has certain rights in the invention. 

This spplication claims the benefit under 35 U.S.C. § 120 of prior U.S. Provisional Patent 
AppKcation Serial No. 60/309,572 filed August 2, 2001, entitted "A METHOD FOR 
10 FREQUENCY ENCODED SPATIAL FILTERING TO ENHANCE IMAGING QUALITY OF 
SCATTERING MEDL«L" 

FIELD OF THE INVENTION 
This invention relates to the field of linear equations, and more particularly to using 
correction filters to enhance solutions to a system of linear equations such as tiie type of 
1 5 equations used in imaging of a scattering medium. 

BACKGROUND 

Imaging of a scattering medium relates generally to a modality for generating an image of 
the spatial distribution of properties (such as the absoiption or scattering coefiBcients) inside a 
scattering medium through the introduction of energy into the medium and the detection of the 

20 scattered energy emerging firom the medium. Systems and methods of this type are in contrast to 
projection imaging systems, such as x-ray. X-ray systems, for example, measure and image the 
attenuation or absorption of energy traveling a straight line path between the x-ray energy source 
and a detector, and not scattered energy. Whether energy is primarily highly scattered or 
primarily travels a straight line path is a function of the wavelength of the energy and medium it 

25 is traveling through. 
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Imaging based on scattering techniques permits the use of new energy wavelengths for 
imaging features of the human body, earth strata, atmosphere and the like ibst can not be imaged 
using projection techniques and wavelengths. For example, x-ray projection techniques may be 
adept at imaging bone structure and other dense objects, but are relatively ineffective at 
S distinguishing and imaging blood oxygenation levels. This is because the absorption coefficient 
of blood does not vary significantly with blood oxygenation, at x-ray wavelengths. However, 
infiared energy can identify the spatial variations in blood volume and blood oxygenation levels 
because the absorption coefficient at these wavelengths is a function of hemoglobin states. Other 
structures and functions can be identified by variations or changes in the scattering coeffici^t of 

10 tissue exposed to infiaied mergy, such as muscle tissue during contraction, and nerves during 
activation: These structures could not be imaged by projection techniques because projection 
techniques are not effective in measuring variations in scattering coefficients. These measures, 
obtainable through imaging based on scattering techniques, such as optical tomography, have 
considerable potential value in diagnosing a broad range of disease processes. 

15 A typical system for imaging based on scattered energy measures, includes at least one 

energy source for illuminating the medium and at least one detector for detecting emerging 
energy. The energy source is selected so that it is highly scattering in the medium to be imaged. 
The source directs the energy into the target scattering medium and the detectors on the surface 
of the medium measure the scattered energy as it exits. Based on these measurements, a 

20 reconstructed image of the internal properties of tiie mediimi is generated. 

The reconstruction is typically carried out using "perturbation methods." These methods 
essentially compare the measurements obtained fiom the target scattering medium to a known 
reference scattering medium. The reference medium may be a physical or a fictitious medium 
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which is selected so that it has properties that are as close as possible to those of the medium to 
be imaged. Selecting a reference medium is essentially an initial guess of the properties of the 
target, hi the first step of reconstruction, a "forward model" is used to predict what the detector 
readings would be for a particular source location based on the known internal properties of the 

S reference medium. The forward model is based on the transport equation or its approximation, 
the difRision equation, which describes the propagation of photons througjbi a scattering medium. 
Next, a perturbation formulation of the transport equation is used to relate (1) the diiSerence 
between ihe measured and predicted detector readings fix>m ttie target and reference, 
respectively, to (2) a difference between the unknown and known internal properties of the target 

10 and refCTence, respectively. This relationship is solved for the unknown scattering and 
absorption properties of the target. The final distributions of mtemal properties are then 
displayed or printed as an image. 

Imaging systems and methods based on scattering techniques, such as optical tomography 
systems, provide a means with which to examine and image the internal properties of scattering 

IS media, such as the absorption and diffusion or scattering coefficients. However, the 

aforementioned imaging systems and methods that recover contrast features of dense scattering 
media have thus far produced results having at best modest spatial resolution. Strategies for 
improving image quality are known (e.g., Newton type), but invariably these are computationally 
intensive and can be quite sensitive to initial starting conditions. 

20 Central to the method of image formation in magnetic resonance imaging (MSI) is that 

. there is a one-to-one correspondence between the frequency of the measured induced current and 
the spatial orientation of the magnetic field gradient. Because the spatial orientation of the 
magnetic field gradient is known, this correspondence pemiits a direct assignment of a measured 
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response to the origm of the signal in space. In effect, the physics of the magnetic resonance 
phenomenon encodes a frequency signature into the measured data that has a known spatial 
relationship with the target medium. More generally speaking, methods of tibls type are known 
as "frequency encoded spatial filtering." 

5 For the foregoing reasons, Hbere is a need for a computationally ef&cient nonlinear 

correction method that is capable of significantly inq>n>ving the quality of solutions to a system 
of linear equations such as reconstructed images of a scattoing medium. 

SUMMARY OF THE INVENTION 
The present invention satisfies this need by providing a method and system for image 

10 reconstruction and image correction that is conq)utationally efficient and improves the quality of 
reconstructed images of a scattering medium. 

In one embodiment of the system and method of the present mvention, reconstructed 
images of a scattering medium are enhanced by: (1) subdividing a target medium into a plurality 
of volume elements and assigning a modulation frequency to at least one of the volume 

1 5 elements' optical coefScients; (2) directing energy into the target medium fix>m at least one 

source during a period of time, and measuring energy emerging from the target medium through 
at least one detector; and (3) processing the measured eaergy emerging firom the target medium 
to reconstruct at least one image, detemiining a frequency encoded spatial filter (FESF), and 
applying the FESF to at least one reconstructed image. 

20 In another embodiment of the system and method of the present invention, a solution to a 

system of linear equations is enhanced by: (1) modeling a target medium into a plurality of 
elements and imposing at least one localized fluctuation into the target medium; (2) measuring 
an output resulting fix>m at least one localized fluctuation; and (3) processing the measured 
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output to reconstruct a result, detenninmg a correction filter, and applying the correction filter to 
tiie result 

The above advantages and features are of representative embodiments only, and are 
presented only to assist in understanding the invention. It should be understood that they are not 
S to be considered limitations on the invention as defined by tiie claims, or limitations on 
equivalents to the claims. For instance, some of these advantages may seem mutually 
contradictory, in that they cannot be simultaneously implemented in a single embodiment 
Similarly, some advantages are primarily sqjplicable to one aspect of the invention. Thus, ttiis 
summary of features and advantages should not be considered dispositive in determining 
10 equivalence. Additional features and advantages of the invention will become apparent in the 
following description, firom the drawings, and fiom the claims. 

BRIEF DESCRIFnON OF THE DRAWINGS 
For a better understanding of the invention, together with the various features and 
advantages thereof reference should be made to the following detailed description of the 
1 S preferred embodiments and to the accompanying drawings therein: 

FIG. 1 is a schematic of an optical tomography system used in accordance with the 
present invention; 

FIG. 2 is an illustrative flowchart describing the method of the present invention; 

FIG. 3 illustrates a target medium; 
20 FIG. 4A illustrates the global spatial correlation between each amplitude map and the 

known spatial distribution of the corresponding frequency in the target mediiun, plotted as a 
function of the modulation frequency for the weight-transform singular-value 
decomposition (SVDWT) algorithm; 
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FIG. 4B illustrates each axiq>litude map's center of mass plotted as a fimction of fm for 
the SVDWT algorithm; 

FIG. S A illustrates the global spatial coireladon between each amplitude map and the 
known spatial distribution of the corresponding frequency in the target medium plotted as a 
S function of fm for the combined SVDWT with an additional matrix preconditioning operation 
(SVDWTWRS) algorithm; 

FIG. SB illustrates each amplitude map's center of mass plotted as a function of fm for 
the SVDWTWRS algorithm. 

DETAILED DESCRIPTION 

10 1. Introduction 

The system and method of the present invention will be discussed in accordance with its 
^plication to the field of optical tomography. It is noted, however, that this methodology 
spplies to a broad range of probl^ns dealing with linear applications in which linear perturbation 
theory is applied to foster a solution such as economics, quality-control, epidemiology, 

1 S meteorology, or the like. 

Image reconstruction methods employ computation-intensive algorithms, which are 
modifications of a standard linear pertuibation approach to image recovery. One of the fiictors 
that has made the development of these algorithms difBcult in the past has been the absence of a 
way to quantitatively characterize the information spread function (ISF) associated with a given . 

20 image reconstruction method. The term ISF used herein refers to the precise manner in which 
the optical coefficients that actually are present in a given location of a target medium are 
mapped into the spatial domain of the image. 
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In flie absence of infoimation regarding the ISF, there is no apparent way of 
systematizing the process of modifying a reconstruction algorithm in response to the observed 
quaUty of its perfonnance. In order to characterize the ISF for a given combination of 
reconstruction algorithm and reference medium the present invention utilizes the techniques 
5 found in magnetic resonance imaging (MRI), which encode a frequency response into 

measurement data that has a known spatial relationship with a target medium. Where the present 
invention differs from MRI is that rather than directly applying this strategy for image fonnation, 
the present invention instead lilies this concept to derive a frequmcy encoded spatial filter 
(FESF) that is then applied to improve the spatial convolution of images previously recorded 

10 using other methods. 

This is accomplished by recognizing that FESFs can be derived by examination of the 
position-dependent temporal frequency spectra obtained from a time series of images whose 
optical properties in each element were assigned dififerent time-varying properties. In the case of 
a perfect imaging method, analysis of the time series would exactly recover the temporal 

IS behavior in every pixel. In practice, spatial convolution is present, in which case the location 
and ampUtude of ttie convolving contrast feature can be determined from exanodnation of flie 
frequency spectrum of the pixel data. However, by assigning temporal properties that are 
uniquely distinguishable among all pixels, precise assignment of image contrast from any one 
pixel to any other is possible. The resulting information is then used as a linear operator that 

20 serves to rearrange (i.e., deconvolve) the contrast features of a recovered image from a test 
medium, thereby improving image quality. Implicit in this scheme is the assumption that the 
spatial convolution defined by the FESF is similar to the convolution present in the image of die 
test medium. In principle, any number of FESFs can be derived and applied as needed. 
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For illustration pmposes, the present system and method is described in further detail 
below with respect to an optical tomography system used to generate images of a targ^ 
scattering medium. Howev^, it will be qjpredated by those sldlled in the ait fbst flie 
methodology of the present invention is applicable in image reconstruction from measured data 

5 based on any energy source (e.g., electromagaetic, acoustic, etc.)* any scattmng medium (e.g., 
body tissues, oceans, foggy atmospheres, geological strata, and various materials, etc.), any 
source condition (e.g., time-independent, time-haimonic, time-resolved) and any physical 
imaging domain (e.g., cross-sectional, volumetric). Accordingly, this methodology can be 
extended to allow for new imaging approaches in a broad range of applications, including 

10 nondestructive testing, geophysical imaging, medical imaging, and surveillance technologies. 
2« . Optical Tomography System 

There are many known imaging systems for collecting the measured data used in image 
reconstruction in scattering media. A schematic illustration of an optical tomography system is 
shown in FIG. 1. This system includes a computer 102, energy sources 104 and 106, a fiber 
15 switcher 108, an imaging head 110, detectors 112, a data acquisition board 114, source fibers 
120 and detector fibers 122. 

The energy sources 104 and 106 provide optical energy, directed through abeam splitter 
1 18, to the fiber switcher 108 and then to each of the plurality of source fibers 120 one at a time 
in series. The source fibers 120 are arranged around an imaging head 1 10 so that the energy is . 
20 directed into the target medium 1 16 at a plurality of source locations around the target 

The energy leaves the source fiber 120 at the imaging head 1 10 and enters the target 
medium 1 16 centered in the imaging head 110. The energy is scattered as it propagates through 
the target medium, emerging firom the target medium at a plurality of locations. The emerging 



wo 03/012736 



PCT/US02/24520 



energy is coUected by flie detector fibers 122 airanged around ihe imaging head 1 10. The 
detected energy then travels through the detector fibers 122 to detectors 1 12 having energy 
measuring devices that generate a signal corresponding to ttie measurement The data 
acquisition board 1 14 receives the measurement signal for delivery to the computer 102. 
5 This process is repeated for each source position so fliat a vector of measures are obtained 

for all of the detectors and source locations. The computer 102 or other suitable processing 
device or hardware is used to process the collected data and reconstruct the image as described in 
detail by tiie methods below. 
3. Method 

10 Figure 2 is an illustrative flow chart describing tibie method of the present inventioiL The 

first step, in accordance with the present invention is to subdivide a scattering medium for which 
the filter fimction will be computed into N small area or volume elements (step 200). Next, a 
sinusoidal temporal variation is assigned to an optical parameter (e.g., absorption and/or 
scattering coefScients) in each area/volume element, with a different firequency to each location 

15 (fu fly • - •> /a) (step 210). The oscillation firequencies (i.e., modulation firequency) in step 210 
are chosen in such a way that every ratio of firequencies /m//n, n /n, is an irrational niraiber. 

The subsequent step involves computing a time series of forward problem solutions: 
lUJk), where j = 1,2,...^ is the detector index and ^ = 1,2,...,*: is the time-step index for the 
resulting dynamic medium with N > 2max(/iO/min(A/m) (step 220). In order to prevent 

20 frequency aliasing in step 220, the time interval between successive states of the medium must 
be small relative to the reciprocal of the highest frequency in the medium. Further, to ensure that 
there will be sufSciently high firequency resolution in the computed time series, the total duration 
of the measurement must be long relative to the reciprocal of the smallest difference between any 
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two assigned fi^encies. The data obtained in step 220 constitute aJx iC matrix of detector 
readings. 

The next step involves solving K inverse problems (i.e., reconstructing the time series of 
tomographic images), one for each set of detector readings computed in stq> 220 (step 230). In 

5 step 230, each column of the Jx matrix of detector readings, in turn, is used to generate the 
left-hand side of the equation 51 = W5x , and the corresponding x is calculated. The data 
obtained in this step constitute an i^x matrix of reconstructed optical parameters - the n n 
row is the time series for the optical parameter in the pixel or voxel. 

Once the image time series is complete, the temporal discrete Fourier transform (DFT) 

10 for each pixel of tiie tomographic images is conQ)iited (step 240). Subsequently, a spatial map of 
the DFT amplitude at each modulated frequency is created (stq) 250). The FESF is determined 
by concatenating the DFTs computed in step 240 into an array or matrix, wherein each row 
corresponds to the DFT amplitude in one image pixel and each column coiresponds to the DFT 
amplitude in one image pixel and each column corresponds to the spatial map of the amplitude at 

15 a particular frequency (step 260). The result, in step 260, is a single linear system, e.g., fia* = 
F^a> where Ha and //a are the iV^x 1 vectors of reconstructed and true absorption coefficients, 
respectively, and F is an Nx matrix that is determined, as described subsequentiy, by 
comparing the matrix of DFT amplitudes conq)uted from the image time series to ttie known 
ideal DFT ampUtude matrix. Deteraiination of = F/ia, i.e., FESF, is accomplished via a 

20 straigfhtforward LU decomposition (i.e., Gaussian elimination). It is noted that a singular-value 
decomposition (SVD) may also be used. Application of the FESF to each reconstructed image of 
the time series is a matter of performing a simple back-substitution (i.e., a spatial deconvolution 
correction to the reconstructed images) (step 270). 
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3 J Forward Model 

The following discussion regarding the Forward Model (i.e., ttie forward problem) is 
provided to elucidate the first step of reconstruction, which is used to predict what fiie detector 
readings would be for a particular source location based on the known int^nal properties of a 
S reference medium. 

As discussed above, typical reconstruction techniques are based on perturbation methods 
that essmtially relate the difiCTence between predicted detector measurements fiom a reference 
medium and detector measures from tiiie target, to solve for the difference between unknown 
properties of the target and known properties of the reference. Accordingly, one of the first steps 

10 in reconstmctipn is to select a reference medium and predict the detector readings by modeling 
or physical measme. Modeling the energy propagation in the scattering medium is done iising 
the transport equation or its q>proximation, the diffiision equatiorL The equations describe the 
propagation of photons through a scattering medium. For a domain having a boundary dA, tins 
is represented by the expression: 

15 V.[/>(r)V«(r)]-A(r)«i(r) = ^(r-r,), reA, (1) 

where u(r) is the photon intensity at position r, ts is the position of a DC point source, and D{r) 
and /^r) are the position-dependent difiEusion and absorption coefBcients, respectively. Here 
the difiusion coefficient is defined as jD(r) = l/{3[jUair) + (r)}, where (r) sflie reduced 
scattering coefficient. Using this equation, the energy emerging &om the reference medium at 

20 each detector location for each source location is predicted. The transport or diffusion equations 
are also the basis for formulating the perturbation or inverse formulation used in reconstruction. 
3.2 The Inverse Formulation 

11 
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The following discussion tegaiding the Inverse Formulation (i.e., the inverse problem), is 
provided to elucidate the second step of reconstructing the time secies of tomographic images. 

As discussed above, reconstruction of a cross-sectional image of the absorption and/or 
scattering properties of the target medium is based on the solution of a perturbation or inverse 
S formulation of the radiation tran^rt or difiusion equation. The perturbation method assumes 
that the composition of the xmknown target medium deviates only by a small amount from a 
known reference medium. Ibis reduces a highly non-linear problem to one that is linear with 
respect to the difference in absorption and scattering properties between the target medium under 
investigation and the reference medium. The resulting optical inverse or perturbation 
10 formulation is based on the normalized difference method and has flie following form: 

W<»^> -5^, +W;^> -SD^Sji,, (2) 

whesrt bfia and 5D are the vectors of cross-sectional differences between flie optical properties 
(absorption and difiusion coefficients, respectively) of a target (measured) medium and of a 
reference (computed or measured) medium used to generate the initial guess; Wj*^ and wj"* are 
15 the weight matrices describing the influence that localized perturbations in the absorption and 
difiusion coefficients, respectively, of the selected reference medium have on the surface 
detectors; and 8ur represents a normalized difference between two sets of detector readings, 
which is defined by the equation: 

20 (8u,),=fil%Ml(„,),. ,=1.2.....^. (3) 



Here, Uris the computed detector readings corresponding to the selected reference medium, Ui 
and Ui represent two sets of measured data (e.g., background vs. target, time-averaged mean vs. a 
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specific time point, etc.) and M is the nimiber of source-detector pairs in the set of 
measurements. 

3.3 Weight Matrix Scaling 

The following discussion regarding Weight Matrix Scaling is provided to elucidate the 
scaling of the weight matrices arrived at in the inverse problem. 

Hie effect of scaling the weight matrix is to make it more unifoim, which can oAea serve 
to improve its conditioning. A scaling qjproach tiiat scales each column of W/^' and Wj^^ to the 
average value of the column vector is used. However, it should be understood that any of the 
known scaling approaches could be adopted. The form of the resulting new weight matrices is: 
v5^»)=ViK*>.R(*). (4) 

where kcmh&fia or A and R^*^ is the normalizing matrix whose entries are: 
1 



. M J = i' 

il(w;*0^ .•.y=i.2.....Ar. (5) 



0 J^h 

in which N is the number of elements used in discretizing the domain A. The resulting system 
equation is: 

Wr^'^ • 6pa + w;°^ • 65 = (6) 

where^ji. =[r^^J* • and 8D=[r^''^ JUd . Note that R^*^ is a diagonal matrix (Eq. 5) 

all of whom main diagonal elemoits are non-zero (Eq. 5); consequently it has a well-defined 
inverse, the computation of which is a trivial matter. 
4. Frequency Encoded Spatial Filter (FESF) 

The following discussion regarding FESF is provided to elucidate determination of the 
FESF, which is used to reconstruct coirected images. 

13 
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The anq)litude "spatial maps" produced by image reconstructioii and tiie computations of 
the DFTs in actuality are strings of numbm, each being the amplitude, at one particular 
fiequency, assigned by the reconstruction algorithm to one of the FEM nodes (i.e., tiie number or 
vertices, or points where three or more elements come together). The entire set of amplitude 
maps can be concatenated into a matrix: 



A,= 



Ai Aa 



An, 

An, 



where N/isihe number of frequencies number of finite elements) and Nn is the number of 
FEM nodes. This is important in understanding the determination of the FESF, because in 
practice the number of nodes invariably is smaller tiian tiie number of elements. Then At (i for 
1 0 "Cages'*) is not a square matrix, but has approximately half as many rows as columns. 

As such, a second matrix can be written, which tells us exactiy where each modulation 
firequoicy actually was pr^ent in the medium: 



B - 



^21 ^72 



\N, 



^2N, 



The matrix Bm must be sparse (i.e., most of its elements are zeroes), because each fm is assigned 
15 to only one of the medium's finite elements. In fact, every row of Bw, which contains hundreds 
or thousands of elements altogether (i.e., Nn = 10^ - 10^, has exactiy three (in the case of two- 
dimensional media) or four (three-dimensional media) elements that are not zero. The niunber is 
three or foiu: due to the use of triangular or tetrahedral elements, so each element is bounded by 
three or four nodes. 

14 
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Ifflie previously described reconstniction process we^ This ideal 

result, however, is not achieved in practice. Thus one must make some type of assumption 
regarding the nature ofthe function that transforms^ Theonemadeintfaepresoit 
invention is that the frequency spectrum present at any one node (i.e., pixel) in the images is a 
5 linear function of the frequencies present at all nodes in the medium. Mathematically, tiiis is 
stated by: TA/ = Bm, where T is a Nn^Nn (i.e., square) matrix. In practice, one wants the 
transformation to go in the other direction, that is, starting fit>m Bm, produce something that is as 
close as possible to die true A;. Thus, computation of the filter that will actually be used in 
practice is accomplished by solving the matrix equation A4 = FBm, where F also is a Nn^Nn 

10 matrix. T and Fare inverses of each other. 

It is noted that the total number of elements in F is smaller than the number in Bm, 
because Nn < Nf. This means that perfect correction will not occur whm applying this method, 
because there aren't enough correction teems to go aroimd. This occurs, not because of the 
assumed linear relation between A, and B^^ but becaxise there is an unavoidable loss of 

1 S information associated with mapping the frequencies in Nf elements into the smaller number, Nn, 
of nodes. 

The FESF that is computed in this way has quality-control utility as a way of quanti^dng 
the accuracy of reconstruction algorithms. The FESF may further be used as an image enhancing 
tool if it is employed in conjunction with data obtained from different experimental media from , 
20 the one used to generate the filter. In this scmario, upon reconstruction of a set of images Ii, I2, 
etc., then, to the extent that the filter function is not strongly dependent on the medimn's 
properties, the spatial accuracy of the reconstruction can be improved by computing FIi, FI2, etc. 
5* Demonstration Results 
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The following example is presented to illustrate features and characteristics of the present 
invention, and is provided solely to assist in explanation of a demonstration of the invention and 
is not intended to be construed as limited ttiereto. 

A demonstration of the utility of the FESF is described in the foregoing example. As 
S discussed below the FBSF has been implied to two different image time series, both obtained 
fiom the same sets of detector readings but enq)loying difiTerent varieties of a reconstruction 
algorithm. In principle the reconstruction methods employed should produce identical results 
since there is no self-evident a priori reason for choosing to use one rather than the other. 
However, application of the FESF method indicates that one variety of reconstraction methods 
10 can produce spatially accurate images of perturbations at any location of the modeled medium, 
and ttie other can not do so. Accordingly, the computed ESF for either algorithm affords a way of 
applying a spatial deconvolution correction to a reconstructed image. 

Figure 3 illustrates a regularly-shaped two-dimensional medium (i.e., the target 
medium). As shown in FIG. 3, the medium is a homogeneous disk of 8 cm diameter, with 
15 optical coefficient values of fl=0.06cm"^ , =Dcm^ For more convenient solution of the 
forward and inverse problems, the mathematical boundary of the disk was extended 0.5 cm 
beyond that of the "physical" medium, as indicated in FIG. 3. The coefficient values in the 
extended region were the same as those of the '^physical" medium. Sixteen equally spaced, unit- 
strength, homogeneous point sources were placed in the mediimi at the indicated positions on the, 
20 physical boundary. 

The numbers of finite elements and nbdes in the indicated mesh are 1604 and 850, 
respectively, and the smallest and largest element areas are 0.026 cm^ and 0.073 cm^ (mean ^ 
standard deviation = 0.040 ± 0.006 cm^). Sinusoidal modulation was imposed on the absorption 
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coefficient in each element. A unique fm was assigned to each, while tibie anaplitude was 
everywhere 0.006 cm"^ (Le., 10% of the mean value). For this preliminary study, Ihe elements* 
scattering coefficients were not moduhited in time. 

To ensure that the resolution bandwidth was smaller than the smallest difference between 
fmS and the Nyquist fiequency was greater than the largest/ii, a time series often thousand sets of 
tomographic detector readings was computed, witii / = 0.01 s. Image reconstruction was 
carried out with two algorithms, both based on an SVD of the image operator matrix. The first 
algorithm used was tihe previously described weight-transform S VDWT method. The second 
reconstruction method - SVDWTWRS - combined SVDWT with an additional matrix 
preconditioning operation, in which each equation was scaled so that all rows of the weight 
matrix had the same sum. 

Two types of anal>^is were performed on the 1,604 DFT amplitude maps produced 
during image reconstruction. First, the global spatial correlation was computed between each 
amplitude map and the known spatial distribution of the corresponding frequency in the target 
medium (ideal result: correlation exactly equal to 1.0 at all frequencies). Second, the coordinates 
of each map's center-of-mass were computed, from which we easily determined its 
displacement from the geometric centroid of the finite element whose a was modulated at the 
corresponding frequency (ideal result: displacement exactly equal to 0.0 at all frequencies). 
These two quantities are plotted, as a function offm (or, equivalently, location in the target 
medium), for the SVDWT algorithm in HGS. 4A and 4B, and for the SVDWTWRS algorithm m 
FIGS. 5A and 5B. The lighter-^lored curve m FIGS. 4A and 4B, and 5A and 5B are derived 
from the unfiltered FT amplitude spatial distributions. The darker cmves are the results obtained 
when the calculations were repeated after we made the best-^)ossible correction consistent with 
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the theoretical model described above, according to which tibe amplitude maps derived torn the 
reconstructed images are a simple linear transfomiation of the true spatial distributions present in 
the target medium. 

Inspection of FIGS. 4A and 4B, and 5 A and SB reveals diat each plotted function exhibits 
5 a qualitative change in bdiavior after the 400^/„. The change is simply a consequence of the 
&ct that the first 400 finite elements all were located in the zone (see FIG. 3) lying between flie 
physical and extended boundaries, i.e., outside the ring of sources and detectors. Closer 
inspection of FIGS. 4A and 4B reveals that both spatial accuracy measures fall particularly far 
fix)m their ideal values for those finite elements corresponding to roughly the 800* througji 

10 1100^/;,. These elements are the ones that lay in the central region of the target medium. That 
is, the SVDWT al^rithm reconstructed images that were strongly distorted spatially, with flie 
absorption coefficient values of the central region significantly displaced toward the surface 
while those of the more peripheral region were recovered with considerably greater accuracy. In 
contrast, the spatial correlation and centroid displacement are considerably more spatially 

IS uniform for the amplitude maps derived from the images reconstructed by the (preconditioned) 
SVDWTWRS algorithm. This is a significant observation, as the two reconstruction variants 
theoretically should yield the same solution when both operate on a given set of detector data. 
Finally, it is seen that in each panel of FIGS. 4A and 4B, SA and SB, most points on the dark 
(corrected images) curve lie closer to the ideal value than those on the light (uncorrected mages), 

20 curve. This demonstrates the possibility that information in the ISF could be used to perform 
post-reconstruction enhancement of the images' spatial accuracy. 

It should be understood tiiat the above description is only representative of illustrative 
embodiments. For the convenience of the reader, the above description has focused on a 
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represeatative sample of possible embodiments) a sample that is illustrative of the principles of 
the present invention. The description has not attempted to exhaustively enumerate all possible 
variations. That alternate embodiments may not have been presented for a specific portion of the 
invention, or that fiirther undescribed alternate embodiments may be available for a portion, is 

S not to be considered a disclaimer of those alternate embodiments. Other appUcations and 
embodiments can be conceived by tiiose witiiout departing fix>m the spirit and scope of the 
present inventiott It is therefore intended, that the invention is not to be limited to the disclosed 
embodiments but is to be defined in accordance with the claims that follow. It can be 
appreciated that many of those imdescribed embodiments are within the scope of the following 

10 claims, and otiiers are equivalent. 
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CLAIMS 

What is claimed is: 

5 

1 . A method of enhancing reconstructed images of a scattering medium, comprising: 
subdividing a first target medium into a plurality of volume elements; 
assigning a modulation frequency to at least one of the volume elmients* optical 

coef&cients; 

10 directing energy into the first target medium fix>m at least one source during a 

period of time; 

measuring energy emerging firom the first target medium through at least one 

detector; 

processing the measured energy emerging fix>m the first target medium to 
15 reconstmct at least one image; 

determining a frequency encoded spatial filter (FESF); and 

applying the FESF to at least one reconstructed image of the first target medium. 

2. The method according to claim U wherein the modulation frequency is assigned 
to an absorption coef&cient of a volume element. 

20 3. The method according to claim 1, wherein the modulation firequency is assigned 

to a scattering coef&cient of a volume element. 

4. The metiiiod according to claim 1, wherein the measured energy emerging bom 
the first target medium is processed by employing a pertuibation metiiod. 

5. The method according to claim 4, wherein the perturbation method employed uses 
25 a forward problem solution to reconstruct the tomogr^hic images. 
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6. The method accoiding to claim 4, wherdn the pertiubation method employed uses 
the inverse problem to reconstruct the tomogrq)hic images. 

7. The mefliod according to claim 1, wherein determining the FESF, comprises: 
computing the temporal discrete Fourier transform of the reconstructed tomogrq)hic 

S images; and 

processing the computed temporal discrete Fourier transform to determine the amplitude 
at a modulation frequency associated with its corresponding volume element 

8. The method according to claim 1» wherein the FESF is applied to at least one 
reconstructed image by performing a simple matrix multiplication. 

10 9. The method according to claim 1, further comprismg: 

directing en^gy into a second target medium from at least one source during a 

period of time; 

measuring energy emerging from the second target medium through at least one 

detector; 

IS processing the measured energy emerging from flie second target medium to 

reconstract at least one image; and 

£^plying the FESF determined from the first target medium to at least one 
reconstructed image of die second target medium. 

10. A system for enhancing reconstructed images of a scattering medium, comprising: 
20 means for subdividing a first target medium into a plurality of volume elements; 

means for assigning a modulation frequency to at least one of the volume 
elements' optical coefficients; 

21 



wo 03/012736 



PCT/US02/24S20 



means for directing energy into the first target medium fi:om at least one source 
during a period of time; 

means for measuring energy emerging fiom the first target medium tiuougjti at 
least one detector; 

S means for processing the measured energy emergmg fix>m the first target medium 

to reconstruct at least one image; 

means for detemiining an FESF; and 

means for applying the FESF to at least one reconstructed image of the first target 

medium. 

10 11. The method accordmg to claim 10, wherein the modulation firequmcy is assigned 

to an absorption coefiScient of a volume element 

12. The method according to claim 1 0, wherein the modulation firequency is assigned 
to a scattering coefBcient of a volume element 

13. The method according to claim 10, wherein the measured energy emerging &om 
IS the first target medium is processed by employing a perturbation method. 

14. The method according to claim 13, wherein the perturbation mettiod employed 
uses a forward problem solution to reconstruct the tomographic images. 

15. The method according to claim 13, wherein the perturbation method employed 
uses the inverse problem to reconstruct the tomographic images. 

20 16. The method according to claim 10, wherein determining the FESF, conq>rises: 

means for computing the temporal discrete Fourier transform of the reconstructed 
tomogrEq>hic images; and 
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means for processing the computed temporal discrete Fourier transform to detemiine the 
amplitude at a modulation frequency associated with its corresponding volume element 

1 7. The method according to claim 10, wherein the FESF is spplied to at least one 
reconstructed image by performing a simple matrix multiplication. 

1 8. The method according to claim 10, further conq)rising: 

means for directmg energy into a second target medium from at least one source 
during a period of time; 

means for measuring energy emerging from the second target medium through at 
least one detector; 

means for processing the measured energy emerging from the second target 
medium to reconstruct at least one image; and 

means for applying the FESF det^mined from the first target medium to at least 
one reconstrocted image of the second target medium. 

19. A program stored on a computer readable medium and executable by a processor, 
comprising: 

instruction code which, when executed by the processor subdivides a first target 
medium into a plurality of volume elements; 

instruction code which, when executed by the processor assigns a modulation 
frequency to at least one of the volume elements' optical coefficients; 

instruction code which, whm executed by the processor directs energy into the 
first target medium from at least one source during a period of time; 

instruction code which, when executed by the processor measures energy 
emerging from the first target medium through at least one detector; 
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instruction code which, when executed by the processor processes the measured 
energy emerging fiom the first target medium to reconstruct at least one image; 

instruction code which, when executed by the processor determines an FESF; and 
instruction code which, when raecuted by fhe processor applies the FESF to at 
least one reconstructed image of tiie first target medium. 

20. A program, according to claim 19, Anther comprising: 

instruction code which, when executed by the processor directs energy into a 
second target medium from at least one source during a period of time; 

instruction code which, when executed by the processor measures energy 
emerging from the second target medium tfarpugjb at least one detector 

instruction code which, when executed by the processor processes the measured 
en^gy emerging from the second target medium to reconstruct at least one image; and 

instruction code which, when executed by flie processor applies the FESF 
determined from the first target medium to at least one reconstructed image of the second target 
medium. 

21. A method of enhancing the solution to a system of linear equations, comprising: 
modeling a first target medium into a plurality of elements; 

imposing at least one localized fluctuation to the target medium; 
measiuing an output resulting from at least one localized fluctuation; 
processing the measured output to reconstruct a result; 
determining a correction filter, and 
£q)plying the correction filter to the result. 
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